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Abstract. 



New simulations are reported of a single bilayer immersed in a liquid solvent, us- 
ing a simple extension of the model of the Max-Plack group as used in earlier work. 
Fluctuation spectrum vel structure factor is dissected in detail and the role of bulk 
fluctuationis is revealed.. We propose to search for protrusions directly where they are, 
i.e. at the solvent-head boundary. In this context new single-point quantities and new 
two-point correlations are introduced and determined for the solvent-head pairs. Most 
unusal shapes are obtained. 



Introduction. 

Simulations of liquid bilayers have been numerous and helpful in understanding 
their properties [1] [2]. In one particular line of research the structure factor S(q), readily 
determined in the simulation, is examined. For extended bilayers under lateral tension, 
S shows the 1/q 2 divergence as the transverse Fourier vector q approaches zero. It is 
a manifestation of the celebrated Goldstone mode and the capillary-wave undulations 
related to it. This capillary- wave contribution, augmented as it may be by a 1/g 4 term, 
dies quickly with raising q. 

Beyond that divergence, the fluctuation spectrum vel structure factor is not widely 
nor fully understood. There is present, as must be, some contribution from local 
compression-dilatation processes present in all liquids as density fluctuations. Also, 
the polar heads in contact with the (polar) solvent may not make a perfect plane but 
may "protrude" into the solvent. The morphology of such "protrusions" and the dynam- 
ics of their formation is not known and has not been determined in any simulation. But 
otherwise it has been convincigly shown theoretically [3] that the protrusions, if indeed 
they exist, ought to modify the hydrophobic repulsion between adjacent bilayers in a 
lamellar stack. In the stack, bilayers are separated by layers of water - whose thickness 
d defines various regimes of the physics of the lamellae. The hydrophobic repulsion has 
been measured in the laboratory and its exponential decay with d is established. 

For such reasons, search for protrusions is desirable and much needed. We do that 
in a new way. We propose to look directly for protrusions where they are, i.e. at the 
solvent-head boundary. 

But first we examine further the structure factor S(q) defined as the height-height 
correlation function. (Section 1). As we have emphasized long time ago[4,5] it is neces- 
sary to go beyond the asymptotic region of the capillary wave divergence and examine 
a wider range of the Fourier variable q. Then the nearest- neighbor peak appears [4] [6]. 

In the search for protrusions we introduced some novel quantities, which are shown 
in Section 2. But above all, instead of searching for their presence in S{q), we attempted 
to find them "where they are" and that was met with success. 

The model we use for simulations has been introduced by the Max Planck group [7] 
and represents the liquid as a collection of spherical particles interacting with various 
Lennard-Jones 6-12 potentials. The amphiphilic (surfactant) molecules are made of 5 
particles - beads connected by chemical unbreakable bonds into a single chain. One 
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terminal bead is the hydrophilic head, the remaining 4 beads form the hydrophobic tail. 
This model also goes by the name of "coarse-grained" model and is a simplification of 
the so-called " Martini force field" [8] . Therefore we use molecular units related to the 
model and the potential used; in particular the collision diameter a and the depth of 
the potential minimum e. Thus the reduced temperature is ksT/e. 
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Section 1. Structure Factor and Bulk Correlations. 

The fluctuation spectrum often termed the "structure factor" S(q) is an average 
height-height corrrelation; the very concept originates from a theoretical picture. In the 
theory, a bilayer or membrane is pictured as an infinitely thin sheet - a mathematical 
surface - embedded in three dimensions, initially flat "at rest" or "at equilibrium". For 
small deviations from that state, the state of the surface is fully given by a function 
h(x, j/). The height h is measured along the z— axis and z = zq ,i.e. h(x,y) = zo is 
the reference flat state. The requirement is that for every x, y the height h(x, y) be 
unique, i.e. h(x,y) cannot be a multivalued function. Fluctuations, due to thermal 
motion, are then described by the correlation /i(Ri)/i(R 2 ) . In a translationally invariant 
system, the latter reduces to a function of the relative distance vector R i2 and that is 
Fourier-transformed to become S(q) on averaging. Here R stands for position vector in 
the (x,y) plane and q = (q x ,Qy)- 

Of course no membrane nor bilayer is infinitely thin; wobbling and undulations of 
the surface position h(x, y) are not the only source of fluctuations. A fuller description 
of the fluctuations in any liquid is provided by the density-density correlation function 

ff(l,2) = <p(l)p(2)>, (1) 

with density p(r), abbreviating r = (x\, yi, z\) to (1) etc. For nearly planar bilay- 
ers, exploiting translation invariance and taking the Fourier transform w.r.to the R12 
variable results in a function of three variables H(zi, z 2 \ q)- This function in inhomo- 
geneous systems with planar interfaces or wetting layers has been extensively investi- 
gated[9,10,ll,12]. 

The height fluctuation spectrum S(q) = (h q h- q ) is but a projection of the density- 
density correlation function: 

p L p L 

S(q) = / dz\ \ dz 2 {z 1 - z )H(z 1 , z 2 ; q)(z 2 - z ) (2) 
Jo Jo 

Besides providing a check on the direct computation, this route is instructive because it 
demonstrates that S(q) is a projection of H which, as we know, contains bulk density 
fluctuations and short range correlations. Indeed in a homogeneous system H reduces 
to the Fourier transform of the radial distribution function g(r) (plus an unimpor- 
tant term). It was obvious[4,5] that, besides the small-g asymptotics of bending-cum- 
capillary waves, the traditional S(q) must contain contributions from ordinary density 
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fluctuations present in all liquids. That it is so the nearest-neighbor peak[4,6] in S(q) 
tells us. But to see it, a sufficiently wide range of q must be included in the computation 



In simulations one deals with particles for which the microscopic density is defined 
as the sum of Dirac delta-functions 



where j is the particle index. This is converted to p obtained as the number of particles 
in a given box (bin) of arbitrary size; formally it is obtained by integration over the 
box size. So far we have limited the sum to particles j G (l,iV) being the beads "a" 
that is the heads of the amphiphilic molecules, also distinguishing the "upper" from the 
"lower" monolayer. In principle one can construct and compute the head-solvent (a-s), 
head-tail (a-t), and tail-tail (t-t) correlations but this has never been done so far. 

It is well known that S(q) diverges at q = + showing in a spectacular way the 
existence of capillary waves as (1/g 2 ), possibly augmented by the bending contribution 
(1/ g 4 ) [13] . What happens with the increase of q is not so well understood. In much of 
work only a small range of q was explored; an exception[4] [6] with larger range revealed 
several features. Fig.l shows recent data, of the same qualitative shape. After the initial 
fall, S goes through a minimum near q 2 ~ 1. and raises to a series of peaks, the first 
one located near q = 2% /a, where a is the collision diameter, our unit of length. Hence 
the first peak is the nearest-neighbor peak, very much like the n.n. peak of the bulk 
structure factor. 

The existence of a minimum presents a temptation to represent S as a sum of a 
decreasing and an increasing term. One would be the asymptotic divergent contribution 
of undulation waves , the other would be the contribution of bulk-like fluctuations, 
culminating in the nearest-neighbor peak. 

An excellent idea is now borrowed from Reference [6] where besides S, the average 
with unity replacing h was extracted from simulation. This was deemed to represent the 
bulk-like fluctuations [6]. Formally, we obtain this function by integrating the density- 
density correlation function H over the ^-variables without the h factors. There results 
the generalized susceptibility, already introduced [12] in the context of undulating liquid 
surfaces. First 
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and then 



X(q)= / dz 1 x(z 1 ,q)= / / dzxdz-. 




2H 



(5) 



Thus 



X(0) = ($^$^exp[zqR im ]) 



(6) 
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In the context of an interface in an external field, x can be interpreted as a susceptibility, 
but here it represents a certain projection of the bulk fluctuations of the liquid bilayer. 

Subtractions, i.a. of S(q) — x(<z)> were also investigated [6]. In Fig. 2 and 3 we show 
the results of such subtraction for the smoother version of S i.e. with S of individual 
monolayers and with floating ("recentered") local reference planes. The susceptibility 
x(q) is fitted to a scaled and shifted Lorentzian about the n.n. peak. The capillary- 
wave asymptote has not been subtracted hence it persists at q — >■ 0. The difference is 
not nil, but even without fine-tuning the shapes show this is the right track - the bulk 
contribution is dominant and x(<?) is a close representation of it. One may think this 
fit remarkable. 

The computation of £ from simulation involves a certain detail [6] about the choice 
of reference plane z = zq from which h(x, y) is measured. For a system bilayer -l-solvent 
confined to a simulation box L x L x L z , the obvious "absolute" coordinate system 
has the corner of the parallaepiped as its origin. Let the position of the surface be 
given by h{x,y). Then the flat reference state is h{x,y) = z ; the issue arises what is 
our choice for zq! [6]. Briefly, it can be the middle (plane) of the bilayer z m id or a 
middle (average) plane of either monolayer, z\ or Z2, each either constant throughout 
the simulation run (i.e. constant in "time") or calculated afresh at each time t (i.e. 
"recentered" or "floating"). Explicitely then 



and similarly for z\ or z^ with the sum limited to one monolayer. In our work we 
used the floating individual reference planes z±, Z2 for each monolayer - in most cases. 
These choices were discussed and investigated in detail [6]. That the choice matters is 
proven by Fig.l which shows S(q) computed in the same run with recentered z\(t) , 
Z2(t), and z m id(t). For the first two S coincide (as might be expected) but the last 
shows a transition from the undulation regime to bulk fluctuation regime much better 




(7) 
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displayed. Such an abrupt change is therefore a universal feature as it was seen with 
two other models of intermolecular forces [6] and now with our version of the "coarse- 
grained" model. The region of q where S raises towards the n.n. peak, corresponds in 
the bulk liquid to distances larger than the n.n. distance and therefore to many peaks 
in the r.d.f. g(r). The three-dimensional Fourier transform g(k) of the bulk liquid g(r) 
is smooth just as S(q) is in that region. 

We have also investigated the shape of a function of q calculated by the same 
recipe as S but for a slab of homogeneous liquid of spherical free particles (a) with 
120000 particles at high density p = 0.9 and with periodic boundary conditions or (b) 
immersed in a larger cube of the same liquid. In (b) boundaries of the slab were ghosts 
so that the slab was an open system with a fluctuating number of particles. For a 
thickness of 6 (and zq constant in the middle of the slab), we obtained S(q) of the same 
shape as the Percus-Yevick 1/(1 — pc(q)) curve shown in Fig. 4. Also S(q) (weighted 
by (h — zq)) and x(q) ( n °t weighted cf. eq.(5-7)) differed very little in shape. See Fig. 5. 
Such structureless background with only nearest-neighbor peak ought to remain if the 
capillary- wave contribution is subtracted from full S(q). This exercise also showed that 
x(q) was very close to these quantities obtained for a homogeneous liquid. 

A technical detail is that when we compare S(q) or any other quantity derived 
from h(x,y) for different reference planes, the difference between two correlations < 
X(q)X*(q) > (where asterisk denotes a complex conjugate) and < Y(q)Y*(q) >, can 
be converted to < (X - Y)(X - Y)* >, but the mixed term < X*Y > + < XY* > 
must be included. 

As mentioned in connection with the theory of scattering off surfaces[14], the av- 
erage < (h(R\) — /i(R,2)) 2 > is invariant with respect to the choice of reference plane 
and is equal to < h(R) 2 > —S(Ri2); the invariance of S(q) follows for q ^ 0. Because 
of imperfect orthogonality of the Fourier coefficients, the S(q) from simulations was not 
invariant. 

There exists a method of computing Fourier coefficients with perfect orthogonality; 
it has been employed [15] for computation of the average area of the fluctuating bilayer 
A > L 2 and in a calculation of the areas and curvatures [16]. An illustrative example 
can be found in the Appendix B. Given N coordinates rj interpreted as Xj,yj, h(xj,yj) 
we choose a number n = 2m + 1 and construct the Fourier finite sum 

h(R) = a + ^ a m cos(q m R) + b m sin(q m R) (8) 
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with the set of twodimensional (/-vectors arbitrarily predetermined, e.g. filling a half 
of a square centered on (0,0). The coefficients a m ,6 m are then determined by the 
least-square fit of h to N > n particle coordinates. The price one pays for the perfect 
orthogonality is the limited number of g's; at n = N the function h(x,y) becomes an 
interpolating polynomial and this is not what one wants. In practice n = 2m + 1 ought 
to be significantly smaller than the number N of points to be fitted. An example of 
computed curvature is given in Appendix B. 

1A. The asymptotics of S(q). 

A quantitative description of S(q) for vanishing q is provided by the capillary-wave 
theory which has been formulated and derived in a variety of ways (see e.g. [10]) and 
predicts S cap (q) = (kT/A)(l/(D+jq 2 ). Here D stand for the external-field contribution, 
A is the nominal (projected) area, (3 = 1/kT, and 7 is the true full macroscopic[10] 
surface tension. On the other hand the very successful mesoscopic bending hamiltionian 
was combined with the above result [13]. In the absence of an external field (D = 0) 

S(q) = (kT/A)(l/( iq 2 + Kq 4 ) (9) 

with k the bending (rigidity) coefficient. 
We rewrite (9) as 

f(z) = l/( 1 z + KZ 2 ) z = q 2 . (10) 

Besides the obvious pole at z\ = there is another pole at Z2 = —j/k. Normally 7 > 0, 
z 2 < and is outside the range of z, z > 0. However, when we extended our simulations 
to the entire bilayer isotherm[4], we found a crossover transition to a floppy state (at 
areas smaller than that of the tensionless state), where 7 is negative. Then the pole Z2 
has moved to the positive part of the real axis and, as a result, S(q) diverges not at 
q = but at q* = ^fz 2 = a/M/k. 

This finding [4] qualifies the standard statement about the divergence of the fluctu- 
ation spectrum vel structure factor S(q). It is only in the extended state with positive 
T that S diverges at q = + . 

Splitting f(z) in partial fractions, we find 

((*! - Z 2 )f( X ) = - (11) 

{z — zi) (z - z 2 ) 

Thus the g -4 dependence is an illusion outside the tensionless state. At the tensionless 
state the two poles merge and the asymptote is a pure 1/z 2 = q~ 4 . 
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The standard expression(9-10), proposed heuristically[13], found a support from a 
theory [17] starting from a microscopic bending hamiltonian. There k is given; 7 appears 
as a derived quantity. That calculation can be applied [2] to the height-height correlation 
S(q) and in all versions invariably the general form (9) results(in [2] as based on[17]). 
Either the intrinsic area of the membrane is assumed constant, or microscopic elasticity 
is added, either there are three areas, four, or only two - always a form (9) is obtained, 
asymptotically at the saddle-point it must be admitted [17, 2, 16]. There seems to be no 
escape from it; unless if the theory would dispose of its basic picture of a single infinitely 
thin mathematical surface. 

The standard expression(9) does not work well as a fitting equation for the simu- 
lation data of S(q) at non- vanishing lateral tension T, even at apparently small q where 
other (bulk-like) contributions are supposed negligible. If V > also g > where we 
denote by g the value of 7 obtained from a least-squares fit to S. Then at q — > + , 7 wins 
and S ~ 1/ ' z = 1/q 2 . But adding a positive Kq 4 = kz 2 term to the denominator makes 
S smaller, i.e. would force the curve downwards with the second derivative negative - 
whereas the experimental S goes gently up above 1/z = 1/q 2 . And k must be positive 
if it is to be a valid bending coefficient. These defects are apparent as elucidated below. 

The apparent difficulties with (9) are displayed in the plot of the inverse, 1/S 
against z = q 2 . The data were for T = 1.35, a = 1.888,L = 40, AT = 1800 and 
importantly V = 1.16. With such large T, 7 in eq.(9) will be positive. Therefore the 
curve following the data points, after starting linearly with positive slope should curve 
upwards because k > and yet it did curve downwards. Then we plotted l/zS(z) 
against z; this is Fig. 6. Here we can identify the very small linear part, 7 + zk, which 
agrees with eq.(9). The point is that this range is very small. Having estimated the 
slope k and the origin 7 we go back and add the line 27 + z 2 k, to the plot of inverse 
1/S; this is Figure 7. The data agree with the theoretical prediction, but within a very 
small range of z = q 2 . 

So the theoretical expression (9) can be used and applied, but its range of validity 
is extremely small - in the Figures essentially up to the first two points. Further points 
are " spoiled" by the non- negligible contribution of bulk-like density fluctuations which 
raise S and put down 1/S. 

In all discussion so far we have used the notion of " the average" . By the symbol < 
. > we mean a thermal average according to a (appropriate) Gibbs ensemble equivalently 
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a "time" average over the simulation run; if t is the index of a successive configuration 
of particles and there are T such steps, then 



At a given step t we can introduce another average, that over particles; an example is 
an average ^-coordinate of all particles 



which we take as the floating reference plane z mi d(t). The notation <> does not refer 
to such subaverages at a given time step. 

The empirical subtractions have not shown any apparent protrusions either in this 
work or earlier [6]. We therefore abandoned any further search for their signature in 
S(q) and examined other possibilites. 
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B. Several Correlations and Protrusions. . 

The density-density correlation H but an example of two-point correlations; we 
can define similar correlation functions for any quantity which can be attributed to a 
position in space or to particle j. Let it be C, its spatial density 

JV 

C(x,y,z) = Y,C(j)5( Xj - x)( yj - y)(z J - z) (14) 
i=i 

Correlations Hcc will share several properties with the density-density correlation, 
notably its symmetry. We have produced several such correlation functions, mostly in 
our search for protrusions. 

But before we move on to these complicated objects (functions of 3 variales) we 
show first results for several one-point functions. 

An important observation is that a protrusion is a local irregularity of the bound- 
ary or transition region between the heads of a monolayer and the adjacent solvent.lt 
has not much to do with the spatial distributon of heads, for which it is but a detail 
in the wing of p a (z). That is why it does not show, apparently, in the shape of S(q). 
Fig. 8 shows an example of the densities in this boundary region, for a low-temperature 
tensionless bilayer. Th solvent density p s (z) falls from its high value of 0.89204 inter- 
secting the left wing of the (Gaussian with very slight asymmetry) density of heads in 
the "lower" monolayer. This merging of heads and solvent is a characteristic feature of 
the tensionless states; it is not always possible to point out to a sharp "boundary"; we 
deal with smooth change. 

The derivatives dp a (z)/dz — p' a (z) and p' s (z) are more specific; maximum or min- 
imum may serve as the defined position of the head-solvent boundary. However, these 
quantities are averages over the plane (x, y) and thus will not uncover any protrusion. 
Better candidates are their fluctuations; off-hand one might suppose that with many 
protrusions forming and vanishing the fluctuation of p' a ought to be larger and con- 
versely. The second moment can be computed easily as an average over a simulation 
run for each z. We do not show these results because the protrusions did not show 
visibly enough. 

The densitites are computed by counting the numbers of given species of parti- 
cles in boxes (bins) thus producing a histogram, but the derivatives are not produced 
by numerical differentiation but by exploiting the exact relation between equilibrium 
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averages, known as the First Yvon Equation. It relates the spatial gradient of the 
first-order distribution function (density) to the average force on that point, e.g. for the 
x— component (d/dx)p(x, y, z) =< F x (x, y, z) >. Here F denotes the force acting on the 
volume element dxdydz about x,y,z. For our purposes this definition is transcribed to 
particle notation and in the simulation run we can make use of the forces, all calculated 
at each time step. Thus 



where F z (j) is the z— component of the total force acting on particle j and the brackets 
denote the ensemble average - in our case, average over time. The Dirac-delta localizes 
the particle j at z and this is routinely averaged over bins of chosen size (0.025 in 
all our runs). Checks with numerical differentiation not only showed agreement but 
also a superiority of the force calculation. The fluctuation of dp/dz is then found from 



The large fluctuation in p' was an interesting and new quantity, but disappoint- 
ing as a tool for detecting protrusions. Eq.(15), the first Yvon equation, suggests an 
examination of the instantaneous and/or partly averaged force and this thread was 
developped; some examples are shown below. 

We have investigated several single-point quantities and their histograms P(C)dC, 
e.g. C = cos(6>) where is the angle of F, with the z-axis, F z = F cos(6>); either as an 
average over all pairs of chosen species, e.g. head-solvent," a-s" , or as the angle of the 
total force on particle "j". For angles, the position vectors are equivalently used. The 
wiggling of the boundary between the monolayer heads and the adjacent solvent ought 
to affect the histogram, along with the most-probable angle and the average angle. 

Another possibility is to search for the coordinates x, y, z of this boundary. On the 
average, we obtain the density profiles and we can take the z-coordinate of the crossing 
of falling p s {z) with raising p a {z) as the definition of the position of the boundary. But 
the wiggling is gone on averaging. 

We have tried more sucessfully to use the midpoint of the a-s position vector of 
a-s pairs. Thus for each such pair x c = (xj+x m )/2 with j e a, m e s or conversely - with 
identical definitions for y c and z c . These c-points may be used to construct a surface, just 
as the positions of the heads are traditionally used to construct the surface representing 
the monolayer. But the distribution of their occurence is also of interest. The histograms 
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P(z c )dz c are shown in Fig. 9 for a low temperature and with/without external field 
which smooths the boundary and ought to damp dynamically the protrusions. The 
external field (linear, derived from a double- well parabolic potential) sharpens P but 
also introduces a hump on top of the expected slight asymetry. 

Having a set of positions x c , y c , z c we can treat agan the ^-coordinate as a height 
and define the new S c as 

S c (q) =< h c (q)h c (-q) > (16) 

with 

h c (qx, q y ) = ^2(z c (j) - z Q ) exp[i(q x x c (j) + q y y c (j)} (17) 
j 

and the sum over all c-points. 

The resulting plot of S c vs q has an entirely new shape. Fig. 10 shows two cases: 
with and without external field. S c picks up some capillary-wave contribution, which is 
correct and this is suppressed by the external field. Otherwise the field appears to have 
little effect on higher modes. Here the nearest neighbor peak does not appear at all. S c 
always shows a definite hump in the interesting region of q G (1, 10) where protrusions 
are expected to enhance the fluctuation spectrun. 

The small-gradient theory predicts [18,2] that the area increment is directly related 
to S, AA = A — L 2 = J2 q A q and A q = S(q) * q 2 . Accordingly, we computed S c (q)q 2 ; 
the greatest increment of area precisely falls in the region of q where protrusions are 
expected, as shown in Fig. 12. For comparison, Fig. 11 shows the same data before 
multiplication by q 2 , i.e. it shows the series of S c (q). The series range from the floppy 
state to the most extended state of the bilayer isotherm, all at T = 1.1. Interestingly, the 
hump in Fig. 12 does not change much with the stretching/compressing of the bilayer; 
the three states in the floppy region visibly differ from the rest by an overall raise at all 

q- 

The idea of protruding heads leaving temporarily their environement and venturing 
out into the neighboring solvent is in fact closely related to the concept of the roughness 
of a surface. Theory of equlibrium roughness (see e.g. [14]) of a surface, in particular 
of a solid, is related to the picture of two points on that surface and on the properties 
when their distance greatly increases. This view is corroborated by our results when 
the correlation S c gave such promising results while the single-point fluctuations were 
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not very productive. 
Summary. 

New simulations are reported with an improved "coarse-grained" model and force 
field. A discussion is given of the bulk contributions to the ubiquitous fluctuation 
spectrum vel " structure factor" . It is shown that in order to give any interpretation 
to the spectrum beyond the tiny asymptotic region, the range of Fourier vector must 
include the nearest- neighbour peak at q ~ 2n/a, at the least [4] [6]. The absence of any 
visible signature of the hypothetical protrusions, led to a search for such a signature 
at regions of space where these are supposed to occur, i.e. at the boundary between 
solvent and heads. Among many possibilities which are discussed, the new study of the 
head-solvent correlations appears very promising. The new correlation function S c and 
the corresponding theoretical area increment are reported and a tentative interpretation 
in terms of protrusions is given. 
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APPENDIX A 

In the Appendix we give some details about the model and simulation procedure. 
The model we use was introduced by the Max Planck group [7] as an extension of the 
very first modelling of amphiphilic molecules as dimers in a solvent[18]. The simulated 
liquid system is made of spherical particles, which are either solvent molecules or beads 
connected in linear chains of small length common to all, of unbreakable (chemical) 
bonds, to form amphiphilic molecules. One terminal bead represents the polar bead 
and the remaining 4 the nonpolar hydrocarbon chain. The intermolecular energy is a 
sum of pair interactions which are cut and shifted 6-12 Lennard- Jones functions. 

In general the symmetrical matrices of parameters require 6 energies e a/ 3 and 6 colli- 
sion diameters o a fi- Practically in all simulation work the collision diameters were taken 
the same for all particles. Often the energy parameters were made all equal[7,18,20,21] 
as well. Then, as one author noticed [20] the molecule was not amphiphilic any more. 
Most likely the sterical effects arising from impenetrability of beads, apparently ensured 
for bilayer stability. 

Our model can be viewed as a simplified version of the "Martini force field" [8], 
also used under the designation of a "coarse grained" model [6]. The essential feature 
of an amphiphilic modecule is the permanent chemical bonding of two antagonistic 
groups - the hydrophilic (polar) head and the hydrophobic (non-polar) tail, most often 
a hydrocarbon tail. We denote the solvent particles as "s", the heads as "a" and the 
tail beads as "t". The intramolecular bonds are confined to the small neighbourhood of 
the bond distance by a potential well with infinite barriers. The intermolecular bonds 
are ss,sa,aa,st,at,tt and we take all cr's equal and all e's equal in two groups; while 
£tt = fat = e s t = e, the other group has a stronger force field: e ss = e sa = e aa = AUG*e, 
with AUG > 1. In this way the polar beads "a" and "s" are endowed with a stronger 
attraction. In all work reported here the value AUG = 2. was used. . The cutoff at 
distance r = 2.5a is common to all, except that the repulsion between s-t and a-t pairs 
is obtained by putting the cutoff at r = 2^/ 6 ^ ( and therefore the shift at +e). 

The Molecular Dynamics simulations were done in the canonical ensemble i.e. at 
constant T, V, N with Nose-Hoover thermostat, over 1-7 millions of time steps for sizes 
from 49000 to 121000 particles, 1800-4500 amphiphiles. 
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APPENDIX B. 



Here is given an illustrative example of the curvature calculation with and without 
the small gradient approximation. As explained in the main text, for 2000 heads a 
reasonable number of fitting coefficients would be a half i.e. 1000 which means 500 
(/-vectors which translates to a square of 31 by 31 i.e. for q x = (2tt / L)n x we can have 
n x = 0, ±1, ±2, • • • , ±n max with n max ~ 15. That is a limitation, for values of L used 
in simulations, nearest-neighbor q nn = 2tv is beyond reach of this procedure. Already 
for n max = 4 we have 81 by 81 matrices. After solving for the least-square Fourier 
coefficients, which has to be done in quadruple precision, the function h(x, y) and its 
derivatives are computed as needed for the computation of integrals, like integral (or 
total) curvature 




(B.l) 



where the p's are the first derivatives of h and the mean curvature "H is given by the 
standard expression in the Monge gauge in terms of first and second derivatives. Same 
expression with "H omitted gives the true area. Small gradient approximations were 
reasonably close as might be expected for computations with small n max . Figure 13 
shows the curvatures plotted against the areas for several n max ; the fractal nature of 
the area might have been expected but the (total) mean curvature also follows the same 
trend. Near n max ~ L the period is of the size of the bead but we cannot reach that. 
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FIGURE CAPTIONS 



Caption to Fig.l 

The S(q) =< h.h > correlations for T = 1.35 L x = 61.25 N ~ 4000,p = 0.89, total 
121000 particles, tail length 8. S for all heads shown with plus signs (reference plane 
as recentered average z m id(t)). Otherwise either monolayer (crosses) with individual 
recentered zi,Z2- The const. /q 4 asymptote is shown with line. The remarkable break 
near q 2 = 0.1 is also present in individual monolayer S provided the reference plane is 
z midi constant or recentered. 

Caption to Fig. 2 

The generalized susceptibility x(q) due to bulk-like fluctuations is fitted to a scaled 
and shifted Lorentzian f(q) = f + a/((q — q ) 2 + b 2 ) and subtracted from the full 
S(q) =< h.h >. Even w/o fine-tuning, this fit proves the dominant role of of bulk 
fluctuations in S(q) =< h.h >. Bare S - plus signs, x - crosses, / - line, and the 
difference - boxes. 

Caption to Fig. 3 

The subtraction of generalized susceptibility x(o) from S(q) =< h.h >. Plotted 
against q: S(q) ( floating reference plane for one monolayer) - with plus signs, x ~ crosses, 
fit to x by a Lorentzian - line, the difference - boxes. Clearly the bulk contribution to 
S follows the shape of x- Arbitrary scale, forcing equality of n.n. peaks. Same data as 
Figs.l and 2. 

Caption to Fig. 4 

Correlation x(q) = < exp[z/c.R] > for the bilayer shown with plus signs and for the 
slab immersed in LJ liquid - with boxes, plotted against q. The identical shapes show 
that bulk-like fluctuations produce x m the bilayer. Bulk Percus-Yevick S = 1 + ph = 
1/(1 — pc) is plotted with the line. Density p = 0.9. 
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Caption to Fig. 5 

Fluctuation spectra for a homogeneous slab of N = 120000 L-J particles at T = 1.1 
density p = 0.89. The reference plane is at z = L z /2. x =< PP(q) > - plus signs, s =< 
h.h > - boxes, sf =< h.h > - diamonds, s for all N particles, sf - open system within 
L z /4,3L 2 /4. The corresponding quantities for T = 5 are shown with stars and differ 
surprisingly little from T = 1.1 data. The ordinate is n 2 + n 2 proportional to q 2 . The 
n.n.peak is damped by the logarithmic scale. 

Caption to Fig. 6 

Inverses l/zS(z) plotted against z = q 2 for bilayer (boxes) and two monolayers (plus 
signs and stars) with common reference plane zq = z m id recentered and monolayers with 
individual reference planes recentered (circles). The straight line is a best fit to all, of 
g + kz, at z ->■ 0. Data are for T = 1.35, L x = 40., p = 0.89, N = 1800, a = 1.88889. 

Caption to Fig. 7 

Inverses 1/S(z) plotted against z = q 2 , with labels as in Fig. 5, with the parabola 
g * z + k * z 2 obtained from the straight line of Fig. 5. Extremely small range of validity 
of Eq.(9) is apparent. But valid it is. 

Caption to Fig. 8 

Densities p(z) and their derivatives p'(z) in the lower monolayer, shown as follows: 
solvent - full line falls as error function, heads - full line and shaded vertically. p' s (z) 
line with small boxes, very accurately a gaussian. p' a (z) - full line (black) with positive 
and negative values. At this tensionless state at T=l.l the boundary between heads 
and solvent is fuzzy and broad. 

Caption to Fig. 9 

Histogram of the ^-coordinate of r c (see text) plotted against z for the lower mono- 
layer; solvent is to the left (lower z) and heads to the right, tails largely outside the 
Figure. No external field - plus signs, with a linear external field - stars; with an ex- 
ternal field also supporting an artificial protrusion - circles. For comparison the density 
profile of heads (scaled by an arbitrary factor) is shown with small boxes and shaded 
by vertical lines. 
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Caption to Fig. 10 

The new correlation S c =< h c .h c > for one monolayer with floating local refer- 
ence plane with and without external field, plotted against q 2 . T = 1.1, 1800 am- 
phiphiles,near tensionless state. A weak external field from a double well parabolic 
potential damps capillary waves. See text. 

Caption to Fig. 11 

The S c =< h c .h c > correlation with h c the (Fourier component of) z-coordinate 
of r c plotted against x = q 2 .The series includes floppy states and most extended states 
near the breajing point. For details see Caption to Fig. 12 

Caption to Fig. 12 

The < h c .h c > correlation with h c the (Fourier component of) the ^-coordinate of 
r c (see text) multiplied by q 2 and plotted against z = q 2 .In the small-gradient theory the 
quantity zS(z) is the (/-component of the area increment cL4 g .The successive curves are 
(in this order) for L x = L y equal to 100. (label" 1"), 97., 73. 5 (label" 3"), 64.0 ("4"), 62.5 
("5"), 61.25 ("6") 61.0 ("7") 60. ("8"). and59.("9"). The last three belong to the floppy 
region. An enhancement in the region 1 < q 2 < 10 can be seen, largely independent of 
the area per head. T = 1.35, 121000 particles, 2000 heads in this monolayer. The area 
per head ranges from 1.7826 to 5.035. 

Caption to Fig. 13 (Appendix B) 

For successively increasing numbers of Fourier coefficients (M=12,40,112,364) total 
(mean) curvature and area increment (over L 2 ) computed from the same data on 4000 
heads by the method of Appendix B. Here A A is plotted against total H. There is no 
saturation in sight. Each point in the cluster of data points results from one configu- 
ration of heads, taken succesively from the same simulation run, at T=l. 35, density 0.89 
L = 39. 2000 heads in one monolayer. 



20 



FIG.2 



3500 I — i 1 1 1 1 1 1 r 




-500 >— 1 1 1 1 1 1 1 1 1 

123456789 

q 




Sun Dec 04 22:12:23 2011 



q 



FIG.8 




11 12 13 14 15 16 



z 



1e-06 



0.01 



0.1 



10 



100 



q**2 



Fig.11 



en 

N 



7000 
6000 
5000 
4000 
3000 
2000 L 
1000 - 




0.05 
FIG.6. 







■ 

mm 


i i 
■ 






1 ■ 


■ 

■ 
■ 


1 


v 

2 6 


o 










G 














1 







0.1 



0.15 

z=q**2 



0.2 



0.25 



0.3 




0.1 0.2 0.3 0.4 0.5 

z=q**2 

FIG.7. 



